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Abstract 

The steady state of a vibrated granular gas confined by a movable piston on the top is discussed. 
Particular attention is given to the hydrodynamic boundary conditions to be used when solving the 
inelastic Navier-Stokes equations. The relevance of an exact general condition relating the grain 
fluxes approaching and moving away from each of the walls is emphasized. It is shown how it can be 
used to get a consistent hydrodynamic description of the boundaries. The obtained expressions for 
the fields do not contain any undetermined parameter. Comparison of the theoretical predictions 
with molecular dynamics simulation results is carried out, and a good agreement is observed for low 
density and not too large inelasticity. A practical way of introducing small finite density corrections 
to the dilute limit theory is proposed, to improve the accuracy of the theory. 

PACS numbers: 45.70.-n,47.70.Nd,51.10.+y 
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I. INTRODUCTION 



The behavior of fluidized granular systems resembles in many cases that of ordinary 
molecular fluids. Actually, it is by now well established that generalized Navier-Stokes 
equations describe quite accurately many of the experimental and numerical features of 



m. 



The justification for 



granular flows, especially at low density and small inelasticity 
this fluid-like description, and derivation of theoretical predictions for the transport coeffi- 
cients appearing in the equations, have been intensively studied for some time. Quite often, 
idealized systems of inelastic hard spheres or disks have been considered. For mono-disperse 
models of this kind, the studies carried out include the derivation of thehydrodynamic equa- 
tions to Navier-Stokes order by using kinetic theory methods 



Carlo simulation of (inelastic) kinetic equations 
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[12IJ . and, very recently, linear response theory 

Interactions between grains are inherently inelastic. As a consequence, the kinetic energy 
of isolated granular systems decreases monotonically in time and, in order to keep them 
fluidized, it is necessary to continuously supply energy to them. A prototypical way of doing 
it is by vibrating one of the walls of the container, usually the one at the bottom. Also often, 
the interest focuses on the bulk properties of the system, i.e. on the behavior of the system 
far enough from the walls, where the governing laws are expected to be independent of the 
details of the boundaries. Then, the most appropriate possible way of vibration for this 
purpose is chosen. In many situations of interest, the criterium for this choice is twofold: 
simplicity and avoiding undesired effects such as the induction of propagating waves into the 
system. These goals are formally achieved in the limits of very high frequency and very small 
amplitude; the former as compared with the typical relaxation frequency of the granular fluid 
next to the wall, and the latter with its mean free path. Additional simplifications occur if 
collisions of the grains with the wall are taken as elastic and if the wall is assumed to move 
with a sawtooth profile 15|, Il6 |. 

A kind of idealized walls often used in kinetic theory and particle simulations are the so- 
called thermal walls. By definition, the flux of after-collision particles leavin g it corresponds 
to a Maxwellian flux with the temperature parameter characterizing the wall |l7j. Therefore, 
the shape of the velocity distribution of the particles moving away from the thermal wall is 
independent of the distribution of the ingoing particles. It is evident that thermal walls also 
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provide a mechanism to compensate the energy dissipated in inelastic collisions, once the 
granular gas in the vicinity of the wall tends to have a temperature smaller than the one of 
the wall. It is then not surprising that thermal walls have also been extensively used in the 
literature of fluidized granular gases since a decade ago [l8, Nevertheless, it must be 
stressed that it is not at all evident that thermal walls correspond to any limit of a rapidly 
vibrating plate. Actually, it has been shown that in some cases the stability of systems 
driven by a thermal wall and by a wall vibrating in a sawtooth way can differ [^(J . 

The granular system considered in this paper is confined between two parallel plates 
in presence of a gravity force. The mission of the one at the bottom of the system is to 
fluidize the granular medium, as in the previous studies mentioned above. On the other 
hand, the wall on the top is floating, in the sense that it can move in the vertical direction, 
being supported by the granular fluid below it. As a consequence, the position and motion 
of the upper wall is interrelated with the state of the granular media below it, and the 
boundary conditions to be imposed to the hydrodynamic equations, following from the 
interaction between the gas and the piston, must be determined in a self-consistent way. 
Here, the steady state eventually reached by the granular gas between the two plates will be 
investigated, using the hydrodynamic description provided by the Navier-Stokes equations 
for inelastic hard spheres. Collisions of the grains with the piston will be modeled as hard 
inelastic collisions. This defines in a deterministic way the mechanical interaction between 
the particles and the movable wall. The question addressed afterwards is how to translate it 
into an appropriate boundary condition to be used in the context of hydrodynamics. Here 
this will be done by means of an intermediate stage in which an exact boundary condition 
for the kinetic theory description of the system is formulated. This condition relates the 
distribution functions of particles leaving the wall and approaching it, by expressing the 
conservation of the particles flux at the piston. 

An additional boundary conditions is required to determine the hydrodynamic fields. It 
can be obtained, in equivalent ways, from the energy flux at the vibrating bottom wall or 
from the global balance of energy in the system. In this way, the theoretical prediction is 
completed, and explicit expressions for the fields with no adjustable parameters are derived. 
In order to verify the accuracy of this description, the predictions are compared with molec- 
ular dynamics simulation results. As in other steady states of granular gases, the range of 
applicability of the theory is restricted to values of the restitution coefficient for gas particles 
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collisions close to unity, because of the coupling between inelasticity and gradients. Under 
these conditions, reasonable agreement between theory and simulations is observed in the 
bulk of the system, i.e. outside the kinetic boundary layers next to the walls. This confirms 
the validity of the hydrodynamic description, including the needed boundary conditions, to 
describe vibrated granular gases in quite realistic situations. 

The plan of the paper is as follows. In Sec. DTI a previously derived [21] stationary solution 
of the inelastic Navier-Stokes equations for a vibrated dilute granular system in presence of 
gravity is reviewed. The main results are explicit expressions for the hydrodynamic fields 
of the system having two arbitrary parameters. Additionally, they involve the height of the 
system. The results are particularized in Sec. II III for the granular gas between two plates 
described above. Also, the boundary effects following from the interaction between the 
gas and the piston on the top are formulated as a condition for the gas-piston distribution 
function at contact. In the same section and in Appendix [Bj it is discussed why this is the 
appropriate starting point to derive the hydrodynamic boundary condition and to introduce 
self-consistent approximations, like those used in Sec. HVl An analysis along the same lines 
of the hydrodynamic boundary effects due to the vibrating wall at the bottom is presented 
in Appendix O Also in Sec. IIV| the derived boundary conditions are used to identify the 
arbitrary constants in the hydrodynamic profiles derived in Sec. HTJ There are no adjustable 
parameters in these expressions. 

The comparison of the obtained theoretical predictions with molecular dynamics simu- 
lation results is carried out in Sec. |V] for two-dimensional systems. It includes both the 
detailed description provided by the hydrodynamic profiles and also some global properties, 
like the average position of the piston and the the balance of the total energy of the system. 
A fairly good agreement is observed, especially if some (small) finite density effects are par- 
tially incorporated into the hydrodynamic description, through the equation of state of the 
gas. The papers ends with a short summary and some general comments. 



II. THE GENERAL ONE-DIMENSIONAL SOLUTION 



In this section, some of the results already discussed in ref. 2]J will be shortly reviewed 
and summarized for the sake of completeness. The system considered is a dilute granular gas 
composed of N equal smooth inelastic hard spheres (d = 3) or disks (d = 2) of mass m and 
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diameter a. The position and velocity of grain % will be denoted by r\ and Vi, respectively. 
The inelasticity of collisions between grains is modeled by means of a constant, velocity 
independent, coefficient of normal restitution a, defined in the interval < a < 1. There 
is an external gravitational field acting on the system, so that each particle is submitted 
to a force —mg e z , where go is a positive constant and e z is the unit vector in the positive 
direction of the z axis. 

For steady states with vanishing macroscopic flow and gradients only in the direction of 
the external field, i.e. the z axis, the inelastic hydrodynamic Navier-Stokes equations of this 



system reduce to [21 1 



dp m 
_ = -nmg , (1) 

2 a ( ar M_ T(W = _ 



nd dz \ dz dz r 

Here, n(z) is the local number of particles density, T(z) the local granular temperature, and 
p(z) = n(z)T(z) the pressure. The temperature is defined from the kinetic energy in the 
usual way, but with the Boltzmann constant set equal to unity. Moreover, k is the thermal 
heat conductivity and \i the diffusive heat conductivity, that is peculiar of granular systems. 
More specifically, the generalized Fourier law giving the heat flux q z in the system is 

dT dn 

*> = -«te-l'te- ( 3 ) 

Finally, C (0) (^) is the cooling rate accounting for the energy dissipated in collisions. Upon 
deriving Eq. (pp), use has been made of the local equation p(z) = n(z)T(z), where p is the 
hydrodynamic pressure, valid in the low density limit. 

The expressions of the transport coefficients and the cooling rate appearing in the above 
expressions can be written in the form 

K(a,T) = «*(a)«o(T), (4) 

/i{a,T) = /i*(a)/i Q (T), (5) 

C (0) («,T) = C>)^y, (6) 

with k (T) and rj (T) being the elastic (a = 1) values of the thermal heat conductivity and 
the shear viscosity, 



d{d + 2) 2 T{d/2) (T^ 1/2 



16(d- 1)tt- 



m 



k (t) = i i- 1 °- {d - i \ w 



, (T) .Mp K)1 /v-«, (8) 

87T 2 

respectively, and 

Pol?) = — - — • ( 9 ) 

The dimensionless quantities k*, fi*, and (* are given in appendix [A] They only depend on 
the restitution coefficient a. 

Equation (j2J) shows the physically evident feature that, in this state, hydrodynamic gradi- 
ents are induced by the inelasticity, through the cooling rate. Consequently, a restriction to 
small gradients, as it is the case in the Navier-Stokes approximation used above, also implies 
a limitation on the value of a for which the theory can be expected to apply The interval 
of values of this parameter for which the theory actually provides an accurate description is 
very hard to determine a priori. 

The system is supposed to be confined between two parallel walls located at z = and 
z = L, respectively. The nature of this two walls will be specified and discussed later on. It 
is convenient to introduce a dimensionless length scale £ by 



£ = \/a(a) 

where 



dz — — + 

X(z') mg 



(10) 



X(z)= [o*- x n{z)Y X (11) 

is proportional to the local mean free path, pi = p(z = L) is the pressure of the gas next to 
the wall located at z — L, and 

32(d-l)ir d - 1 (*(a) 



a a = 



(d + 2) 3 r 2 (d/2) [«*(a) - n*(a)] ' 
The £ coordinate is a monotonic decreasing function of z, varying between 



(12) 



C M = £(z = 0) = v^)^- 1 (n z + -g-J (13) 



and 



U = t(z — L) — yfaja)^±. (14) 

In Eq. ( Tl3|) . N z denotes the number of particles in the system per unit of section W (length 
or area) perpendicular to the external field, N z = N/W. It must be noted that the variation 
interval of £ depends on p L and N z , but not on the value of L. 
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The physical meaning of £ can be illustrated by realizing that it is proportional to the 
local pressure at the corresponding height z. This follows from Eq. (JT]), that leads to 



p(z) = mg / dz' n(z') + p L 



mg 



^Ja(a)a d - 1 



(15) 



By substituting Eq. (CQ) into Eq. (j2J) and doing the change of variable defined in Eq. (1101) . 
it is obtained: 

f)2 T l/2 fYpX/2 

£T 1/2 = 0, (16) 

with 



b{a) 



2k* - if 
2(k* - //*) 



The general solution of the above differential equation reads 22] 



where A and B are constants to be identified from the boundary conditions, 

b(a) - 1 ii* 



(17) 



(18) 



via) = 



4(k* 



fJL* 



(19) 



and I v and K v are the modified Bessel functions of first and second kind, respectively 22]. 

Therefore, for the system being considered, the pressure and temperature profiles are 
given by Eqs. (TI5T) and ( 1181) . respectively. Consequently, the density profile is 

P(0 _ mg i 



Finally, the transformation from the £ coordinate to the z one is given by 



(20) 



^Jll{a)a d ~ 



n(£ 



(21) 



that follows from Eq. ([TO]) . It is worth to mention that the presence of the diffusive heat 
conductivity /x in the above expressions is not at all irrelevant. Predictions implie d by its 
existence have been checked both by particle simulations and experimentally 



23 



24 



25 



26J. 



To proceed any further, the boundary conditions of the system at the top (z = L) and the 
bottom [z = 0) must be specified. This will be done in the next section. 
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III. CLOSED SYSTEM WITH A PISTON. THE KINETIC BOUNDARY CONDI- 
TION. 



In order to maintain the system fluidized, it will be assumed that energy is being contin- 
uously supplied to it through the wall located at the bottom, and that this is achieved by 
vibrating it. The simplest possible way of vibration will be considered here, namely with 
a sawtooth velocity profile. This means that all the particles colliding with the wall find 



it with the same upwards velocity Vyy 



16]. Moreover, the amplitude of vibration of 



this wall is taken much smaller than the mean free path of the grains in its vicinity. As a 
consequence, the position of the wall can be taken in practice as fixed at z = with very 
good accuracy. Finally, since the main reason to introduce this vibrating wall is to keep the 
granular matter fluidized, collisions of particles with it will be considered as elastic, for the 
sake of simplicity. Of course, all the above corresponds to a very idealized wall that can not 
be fully implemented in actual experiments. 

Next, the upper boundary condition must be specified. The case of an open system 
(L — ► oo) was studied in ref. 2jJ. Here, a different physical situation will be investigated. 
It will be considered that there is a movable lid or piston on top of the gas, as illustrated in 
Fig. [TJ The piston has a finite mass M, contrary to the vibrating wall at the bottom that 
is taken infinitely massive. The piston can only move in the z direction, remaining always 
perpendicular to it, i.e. parallel to the bottom wall. Its position and velocity will be denoted 
by Z and V z , respectively, so that L corresponds to the average value of Z in the steady 
state. There is no friction between the piston and the lateral walls of the container. 

Collisions of particles with the piston on the top are smooth and inelastic, with a velocity- 
independent coefficient of normal restitution ap, < atp < 1. Therefore, the vector compo- 
nent, v±, of the velocity of the particle perpendicular to the z axis remains unchanged in a 
collision with the piston, 

v'± = v ± , (22) 

where the prime is used here and henceforth to denote after-collision quantities. On the 
other hand, when a particle with a component v z of the velocity collides with the piston 
being V z the velocity of the latter, these values change instantaneously to 

v' z ^v z -^^(l + a P )(v z -V z l (23) 
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FIG. 1: Sketch of the system studied in this paper. 



v = v z+ m 

z z ^ m + M 



l + a P ){v z -V z ). 



(24) 



Therefore, in the collision the total momentum is conserved, the relative velocity g z = v z — V z 
changes to g' z = —apg z , and there is a variation of the total kinetic energy given by 

mM 



AE 



(25) 



2(m + M) 

For ftp < 1 there is a loss of energy in every collision between a particle and the piston. The 
change in the ^-component of the momentum P z of the piston in a collision is 

mM 



AP Z = M{V' Z - v z ) 



1 + a P )g z . 



(26) 



m + M 

Since a collision of a particle with the piston is only possible if v z > V Z) it is AP Z > for 
all the collisions, indicating that momentum is continuously transferred from the gas to the 
piston. A relevant relationship to be used in the following is 



dv' dV' = apdv dV z . 



(27) 



For later use, it is convenient to consider also the so-called restituting collision corresponding 
to the velocities v and V z . It is defined by the velocities v* and V* leading as a consequence 
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of a collision to v and V z . Their expressions are obtained directly by inverting Eqs. fl22l) -(!2 

v* ± = v ± , (28) 

M l + a P 
m + M ap 

V; = V Z + 1 -±^g z . (30) 

m + M a P y K ' 

Also, it is g* = cfp l g z and 

dv*dV z * = a7 P x dvdV z . (31) 

The question now is how to translate the above collision rules into one or more boundary 
conditions for the description of the granular gas below the piston. To discuss this point in 
some detail, let us introduce the two-body distribution function for the piston and the gas, 
&(x, Z, V z , t), defined for an arbitrary state as 

$(cc, Z, V z ,t) = N J dx 2 ... J dx N p(x, x 2 ,..., x N , Z, V z , t), (32) 

where p(xi, x 2 , ■ ■ ■ , x^, Z, Vz, t), with xi = {r^Vj}, is the distribution function for the 
system composed by the piston and the iV grains at time t, and x = {r,v}. Therefore, 
<&(x, Z,V z ,t) is proportional to the probability density of finding the piston at height Z 
with velocity V z and a grain at position r with velocity v, at time t. It is normalized as 

//>oc poo 
dx dZ dV z $(x, Z, V z ,t) = N. (33) 
JO J-oo 

The one-particle distribution function of the gas, f(x,t), can be obtained from $ by 
integration over the piston position and velocity, 

poo poo 

f(x,t)= dZ dV z ${x,Z,V z ,t). (34) 

JO J-oo 

Similarly, the probability distribution for the piston, F(Z, V z ,t), is given by 

F{Z,V z ,t) = ±J dx$(x,Z,V g ,t), (35) 

and it is normalized to unity. No reference to any particular state is involved in the above 
definitions. 

For initial conditions in which all the particles are located below the piston, the two-body 
distribution $ at arbitrary later times can be expressed in the form 

$(x, Z, V z , t) = &(Z - z)$ (x, Z, V z , t), (36) 
10 



where Q(x) is the Heaviside step function defined by Q(x) = 1 for x > and Q(x) = 
for x < 0. The function <$>q(x, Z,V z ,t) is not defined by Eq. fl36|) for 2 > Z, and it can 
be considered as being regular everywhere as well as its derivatives, without restriction. In 
order to discuss the form of $o when the position of the particle is taken next to the piston, 
it is useful to decompose it in the form 

$ {x,Z,V s ,t) = $+(x,Z,V g ,t) + $-(x,Z,V z ,t), (37) 

with 

$+(aj, Z, V g , t) = G(g z )$ (x, Z, V z , t) (38) 

and 

(a;, Z, V z , t) = G(-g z )$ (x, Z, V s , t). (39) 

In the last expression, it is Q(—x) = 1 — Q(x). 

Conservation of the flux of particles at the piston implies that 

$ + (r, v*, Z, V z *, t)dv* dV z * \g* z \5(Z - z) = $_(r, v, Z, V z , t)dv dV z \g z \5{Z - z), (40) 



where v* and V* are the restituting velocities defined by Eqs. (1281) - (l30l) . This is an 
exact relationship, valid for arbitrary density, that can be derived also by starting 
from the evolution equation for <&(x, Z,V z ,t) following from the Liouville equation for 
p(xi, Xo, . . . , xn, Z, Vz,t). Equation (jlDl) is obtained by isolating the singular terms at 



z = Z [27| . When Eqs. fl28l - fl3Tl) are employed into Eq. (jlO]) . it can be reduced to 

$_(r, v, Z, V z , t)5(Z -z) = a P 2 $ + (r, v*, Z, V*, t)5(Z - z). (41) 

Combination of Eqs. fl37|) and (j4ip yields 

$ (x, Z, V z , t)5(Z -z) = {l + otp 2 bp l )®+(x, Z, V z , t)5{Z - z). (42) 

Here bp 1 is an operator acting on the velocities v and V z to its right, replacing them by the 
precollisional values given by Eqs. (l2~8~]) - (l3~0~]) . 

Equation (l4"2j) shows that $0 is fully determined at z = Z if $ + is known at the same 
position. Consequently, upon introducing simplifications or approximations on the value of 
$0 a t the piston, i.e. on $ 5(Z — z), they must refer only to either § + 5(Z — z) or $_5(Z — z). 
Otherwise, the exact relationship given by Eq. (j42H may be violated. A natural question in 
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this context is how relevant is that relation in practice. In other words, is any fundamental 
physical property possibly lost if Eq. (fl2j) is not satisfied in a given approximate description? 
The answer to this question is affirmative as it will be seen in the next section. 



In the limit of a very di 



the Boltzmann equation 271 ]. is to assume that 



ute gas, a simple approximation, similar to the one leading to 



Z, V z , t)S(Z -z) = f(x, t)F(Z, V z , t)Q{g z )5{Z - z), (43) 

therefore neglecting all the correlations between the piston and the particles colliding with 
it before the collision. Using the above approximation into Eq. (j4"2l . it is found that 

$ (x,Z,V z ,t)6(Z-z) = f(x,t)F(Z,V z ,t)5(Z-z) 

+Q(-g z )(a p 2 b p 1 - l)f(x,t)F(Z, V z ,t)5(Z - z). (44) 

Upon deriving this equation use has been made of the relations 

b P l Q{g z ) = Q{gl)b P l = &(-g z )b P \ (45) 

following from the definition of bp 1 and the relationship between g z and g* z . The physical 
meaning of Eq. (|4"4"|) is evident: correlations between the piston and the particles in its 
neighborhood are created by the collisions. 

In this section, only the kinetic boundary condition for the piston at the top of the 
system has been considered. The analysis of the vibrating wall located at the bottom is 
much simpler, and it will be discussed later on. 



IV. HYDRODYNAMIC BOUNDARY CONDITIONS 

In the following, attention will be restricted to the macroscopic steady state described in 
Sees. [II] and IIHI Then, the height of the system L there corresponds to the average position 
of the piston, Z, in the kinetic theory description, while the average velocity of the piston 
is required to vanish, 

i /» /*oo roc 

V z = - / dx / dZ dV z V z $ st (x,Z,V z ) 



N 

POO POO 

dZ dV z V z F st (Z,V z ) = 0, (46) 

JO J -oo 
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where the indexes st are used to refer to properties of the system in the steady state. 
Although the above property can be accomplished in other ways, here the simplifying as- 
sumption that § s t(x, Z, V z ) and, therefore, F st (Z, V z ) are even functions of V z is made. This 
is consistent with the results from molecular dynamics simulations to be reported later on. 
Also, it must be kept in mind that for the state we are considering, $ st depends on the 
position of the particles only through the coordinate z, although it is not made explicit in 
the notation. 

It is interesting and illuminating to compute the force F that the granular gas makes on 
the piston in this state. Because of symmetry reasons, the force only has component, that 
is computed in Appendix [B] with the result 

F z = Wn L T LtZ , (47) 

where Ul is the number density of the granular gas next to the piston, i.e. m — n(z — > L), 
and Tl jZ is a temperature parameter of the gas in the same region defined from the z- 
component of the velocity or, equivalently, proportional to the zz component of the pressure 
tensor, 

If* r /*oo /*oo 

T LtZ = —jdzJdvj dZ J dV z mv 2 z %, st (x, Z, V Z )5(Z - z). (48) 

In the Navier-Stokes approximation used in Sec. [U Tl, z coincides with the temperature 
of the gas at the piston, Tl, and Eq. (T47j) is the expected result, since it agrees with the 
hydrodynamic interpretation of the pressure tensor. In its derivation, a crucial role is played 
by the kinetic boundary condition given in Eq. (j4ip . as it is shown in Appendix [Bl Still 
more, the specific form of $-f, s t is not relevant, as long as be consistently derived 

from it. Otherwise, the pressure of the dilute gas defined as the product of the local density 
times the local temperature, would not agree with the scalar defined from the hydrodynamic 
pressure tensor, characterizing the internal forces in the fluid. 

In order to develop a consistent theory, it is then convenient to express the hydrodynamic 
fields of the granular gas in the vicinity of the piston in terms of only $+, s t, instead of the 
complete distribution $o,st- The expression of the number density introduced above in 
terms of <3>o,st reads 

/p /*oo /"OO 

dz dv dZ dV t $w{x,Z,V z )8(Z - z). (49) 
J JO J-oo 



13 



By means of the decomposition formulated in Eq. (|37|) and using the boundary condition 
(jUJ), the above expression can be put in the form 

n z = i±^n« (50) 

Ctp 

with 

dz dv dZ dV z $ + , st (x,Z,V z )5(Z-z). (51) 

J JO J-oo 

Also, it is easily verified that the local velocity flow vanishes at the piston in the steady 
state, as it should. Finally, consider the temperature of the gas next to the piston, Tl, given 
by 

\n L T L = J dz J dv J™ dZ J°° dV z ^<S> , st (x,Z,V z )6(Z - z). (52) 

It is convenient to distinguish between the perpendicular and z contributions to this tem- 
perature, 

^H^ + ^TU, (53) 
where Tl jZ is defined in Eq. ( H8l) and 

j f f poo poo 

T L ,± = (d _ 1)wj . J dz J dv J dZ J dV 2 mvl$ 0jSt {x,Z,V z )8{Z-z). 



(54) 



The boundary condition (T4T1) leads directly to 

Tl,± = T^L (55) 

with 

I p P POO POO 

T tl = m dz dv dZ dV z mv 2 ± $ +tBt (x,Z,V z )5(Z - z). (56) 

(d — l)n L J J Jo J-oo 

A more involved calculation gives 

_ a P (M - m) (+) (+) 

Ti ' 2 " M — a pin Tl - ~ H * ' (5?) 

where 



j P P POO POO 

T$ = -rrr dz dv dZ dV z mv 2 z $ +>st (x, Z, V Z )5(Z - z), (58) 
n L J J Jo J-oo 

H ( z +) = 2apMm [dz [dv [°° dZ [°° dV z v z V z $ +>st (x, Z, V Z )S(Z - z). (59) 
(M — apm)n L J J Jo J-oo 
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At this point, a simplifying hypothesis on the pre-collisional two-body distribution at contact 
is made. The associated marginal velocity distribution is approximated by a product of 
Gaussian distributions, namely it is assumed that 

POO 

dz / dZ$ + , at (x,Z,V z )5(Z-z) = fit\v)Pst{Vz)Q{g z ), (60) 
Jo 



where 



1/2 

= 24 + WbK) ( ™ ) e ( Gl i 



m 



(d-l)/2 



VMB{V±)= \2^) e_2Ti,± ' (62) 

/ M \ 1/2 i^vf 

P »<^K^>) e " T '' (63) 

Although in the same spirit, Eq. (160]) is a somewhat stronger assumption than the partic- 
ularization of Eq. (1431) for the steady state under consideration. The presence of and 
T^~] in Eq. (16T1) and of X^i in Eq. (1621) is required by consistency with the previous results 
in this section. Moreover, it will be assumed in the following that Tl,± = Tl for the sake 
of simplicity. Because of Eq. (|53|) this implies that also Tl z = Tl and, therefore, all the 
diagonal components of the pressure tensor of the gas in the vicinity of the piston are the 
same, consistently with the Navier-Stokes approximation. 

Using Eqs. (l60l) - (|63l) it is obtained that = and, therefore, Eq. (1571) reduces to 

T l = Tl , 1 = ^^t£. (64) 
M — otpm 

Equations (1501) and (1641) relate the properties of the flux of grains reaching the wall with the 
local hydrodynamic fields of the granular gas next to the wall. They will be employed now 
to derive an expression for the heat flux at the piston, Ql. Using standard kinetic theory 
arguments, this quantity can be computed from the variation of kinetic energy of the grains 
colliding with it, namely it is given by the energy of the flux leaving the piston minus that 
of the flux reaching it. Then, it can be written as 

Q L = j 'dz j ' dv J dZ J dV z ^(v 2 z - v' z 2 )g z <5> + , st (x, Z, V Z )5(Z - z). (65) 

This expression is easily evaluated in the Gaussian approximation given by Eqs. (I601 - (l63l) . 
with the result 

V2 M /i , a\ 1/2 //l/f ™\3/2 r (i + a p )(i + 0)M" 



2 y /z M (l + Qy'* (M-aprn 



-Km I M + m V «p 



(M-apmX 
\ M-m ) 



2(M + m) 



n L Tl'\ 



(66) 
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where 

__ mTp a>p(M — m)mTf 



MT^ +) i M ~ a P m)MT L 



(67) 



The calculation of the power injected into the system through the vibrating wall at the 
bottom is much more direct and it is outlined in Appendix O There, it is shown that the 
heat flux at this wall is given by 

Qo = vwPo, (68) 
with p = p(z = 0) = p(£ = £m) is the pressure of the granular gas in the region just above 



the vibrating wall. Equation (l6"8l was proposed in refs. |15j and 



and used many times 



in the literature since then. To get the value of p , first note that Eq. (fl5|) gives 

Po = Pl + mg N z . (69) 

The value of pl can be obtained by means of Eq. (1471) . after substituting T Lz by T L 
accordingly with Eq. fl64l) . Mechanical equilibrium of the piston in the steady state requires 
that F z = Mg Q and, therefore, 

Pl = n L T L = -yy- ■ (70) 

Once is known, the values of and ^ m defined in Eqs. ( |T3l) and ( JT4l . respectively, can 
be determined. The next step is to identify the constants A and B introduced in Sec. [TTJ 
and present in the expressions of the density and temperature profiles, Eqs. (|T8l) and (1201) . 
To do so, two boundary conditions are needed. They are provided by requiring the values 
of Ql and Qo given above to agree with the hydrodynamic expression of the heat flux q z 
given in Eq. ([3]), particularized for each of the wall boundaries, i.e., 

Q L = lim q z (z), (71) 

z— >Z 

Qo = lim ?*(*), (72) 

z— >0 



with 

dT 

q z (z) = -«o 



K - H J^j -fi mg 



(73) 

3y using Eqs. (ITDjl . (IT51) . and (ITU1) . as well as properties of the modified Bessel functions 



221 ]. the above expression becomes 



q z (0 = 2m 5o ^ («* - /Oe 1- " t^-i(0 " 5/^-i(0] • (74) 
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Making in this equation £ = £ m and equating the result to the expression of Ql in Eq. (J66 
it is obtained: 



where 



e(a, otp) 



B_ _ Iy-i(£ m ) - e(a, ap)I u {£„ 
A Kv-^Zm) + e(a, a P )K v (£, 

I6y/2(d - l)^-^ 2 M /! + . 
d(d + 2) 2 T(d/2) M + m \ ^p~ 



(75) 
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(76) 



2(M + m) 

Next, Eq. (|71jl is particularized for £ = £m and afterwards put equal to the right hand 
side of Eq. ( 1681 . This leads to 



AI (t -i(t M )-BK (l _ 1 (Z M )=£ 



1 /2 



(77) 



d(d + 2) 2 r(d/2) v 4(a)(/- - M *) ' 

An alternative to one of the two boundary conditions used above would be to require 
that the total power dissipated in the system due to the inelasticity of collisions be the same 
as the net heat flux injected in the system through the boundaries. This leads to a relation 
between A and B that is a combination of Eqs. (!75|) and (!77|) . consistently with the fact 
that the condition follows directly from the hydrodynamic equations. This can be realized 
by noting that the energy balance equation in the steady state reads (see Eqs. (J2J) and (j3D) 



dq z 
dz 



whose integration between z = and z = L gives 

d ^ 
Qz,l — q z ,o — — - 



dzp{z)({z). 



(78) 



(79) 



The right hand side of this equation is clearly identified as the energy dissipated in the 
system per unit of time and section due to the inelasticity of collisions. 

To close the identification of the hydrodynamic profiles in the system, a prescription to 
determine the parameter <fr defined in Eq. f l67|) is needed. This in turn calls for an expression 
for the temperature of the piston Tp. It is well known that a peculiar feature of granular 



30, 



31| , that in the present case manifests 



gases is the violation of energy equipartition 
itself by the difference between Tp and Tl- This has been discussed with some detail in {32]]. 
For values of the restitution coefficient a close to unity, a good approximation to the value 
of 4> is obtained by taking T P /T L equal to unity. 
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A relevant quantity characterizing the macroscopic state of the system is the average 

height of the piston, L = Z . The theoretical prediction for it is obtained from Eq. (1211) once 

the values of £m and £ m have been obtained and the density profile is known, 

1 ft" d£ . , 

' (80) 



V. MOLECULAR DYNAMICS SIMULATIONS 

To test the theoretical predictions discussed in the previous sections, molecular dynamics 
(MD) simulations of a two-dimensional (d = 2) system of inelastic hard disks have been 
performed. To avoid undesired boundary effects, periodic boundary conditions with period- 
icity W were used in the direction perpendicular to the z axis. In all the simulations to be 
reported in the following, the grains were initially located in a square lattice with a Gaus- 
sian velocity distribution. The initial position of the piston was slightly above the highest 
layer of grains. Then, the system was evolved in time accordingly with the mechanical rules 
governing the dynamics of the particles and the piston, and it was observed that, for values 
of the restitution coefficient a between 0.9 and 1, a steady state with only gradients in the 
vertical direction and no flow field was reached. Actually, in order for the system to reach 



3e not too large, since otherwise transversal inhomogeneities 



this state, its width W has to 
of the type discussed in refs. 



20| and [33| develop into the system. The results presented 
below have been time averaged once the system is in the steady state, and also over sev- 
eral trajectories. More precisely, once the stationary state was reached, its trajectory was 
followed for 4000 collisions per particle, and 10 different trajectories were used in each case. 

To carry out a systematic analysis of the theoretical predictions, it is necessary to reduce 
the number of parameters of the system by fixing some of them. In the present study, 
the values N = 420 and W = 70a (N z = 6cr _1 ) have been used in all the simulations. 
The reason is that the dependence on N z of the theoretical results follows trivially once 
the hydrodynamic description is assumed to hold, and analysis of the simulation results 
indicates that the above values are appropriate for the purposes here, in the sense that a 
hydrodynamic behavior is observed. Moreover, only results for values of a in the interval 
0.9 < a < 1 will be reported, the reason being that for smaller values of the restitution 
coefficient of the gas, qualitative and quantitative strong deviations from the theoretical 
predictions were found. This is not surprising, since the Navier-Stokes approximation is 

18 



expected to fail beyond the weak dissipation limit, as a consequence of the coupling between 
gradients and inelasticity, as pointed out in Sec. [Til 



A. Hydrodynamic profiles 

In Figs. [2] and [3j the dimensionless pressure (p*) and temperature (T*) profiles, respec- 
tively, are plotted as a function of the length scale £ defined in Eq. ( TIP]) , for two different 
combinations of the parameters M/m, a, and ap, as indicated in the figures themselves. 
Remind that £ is a decreasing function of z. The dimensionless fields are defined by 

T* = ^, = (81) 

It is easily seen that in these units, the theoretical predictions for the hydrodynamic fields 
become independent of the velocity vw of the vibrating wall and also of the acceleration go. 
Of course, this applies as long as the system is fluidized and its density low everywhere. 

Consider first the pressure field. The theoretical prediction is given by Eq. ffl5|) and does 
not involve the parameters A and B appearing in the temperature and density profiles. It 
is represented by the solid lines in Fig. [2J In the simulations, the £ coordinate has been 
measured directly by using the discrete version of Eq. (fIU|) . On the other hand, the pressure 
at each value of z has been obtained from the data for the density and the temperature at 
the same value. In the low density limit considered in the theory, it is p* = n*T*, where 
n* = nav^/go- The results for the pressure represented by empty circles in the figures have 
been computed in this way. A good agreement is observed between theory and simulation in 
the bulk of the system, i. e. outside the boundary layers next to the upper and lower walls. 
The results can be considered as satisfactory, especially taking into account that there are no 
adjustable parameters. Nevertheless, it is true that a systematic although small deviation 
is observed. Equation (IT51) is a quite general result, which only requires for its derivation 
the restriction to the Navier-Stokes order approximation, without any particular expression 
of the transport coefficients or of the (local) equation of state. Thus it seems sensible to 
check whether the origin of the observed discrepancy lies in the way in which the pressure is 
computed from the simulation data, namely by using the ideal gas equation of state. For this 
reason, the local pressure has also been calculated from the MD data by using the equation 



19 




FIG. 2: Dimensionless pressure (p*) profiles for two steady states of the system of inelastic hard 
disks depicted in Fig. [TJ The solid lines are the theoretical predictions given in the main text, while 
the symbols are MD simulation results. The empty circles have been computed using p* = n*T*, 
while the triangles have been calculated by means of Eq. (|82|) . The values of the parameters of 
the system are M = 48m, a = 0.94, ap = 0.99, vyy = for the top figure, and M = 92?n, 

a = 0.98, ap = 0.8, vw = 3y/go<j, for the bottom one 



of state proposed by Grossman et al. (l8| for hard disks at finite density. 



p = nT 7 ^^ , (82) 
n c — n 

where n c = 2a~ 2 / y/3 is the maximum packing number density. The pressure values obtained 
by this procedure are represented by empty triangles in Fig. [2j Now, a much better agree- 
ment is obtained. It is concluded that density corrections are identifiable in the pressure, in 
spite of the density being quite small. Actually, the same values of the pressure are obtained 
if instead of Eq. (1821) . the second virial approximation for the equation of state of a gas of 
hard disks is used. 

The scaled temperature profiles shown in Fig. [3] also exhibit a good agreement between 
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FIG. 3: (Color online) Dimensionless temperature profile T* for the same two steady states as in 
Fig. [2j The solid line is the theoretical prediction given by Eq. (|18j) with the values of A and B 
determined in Sec. IIVI The empty circles are MD results and the dashed line (red online) is the 
best fit to Eq. (|18p by varying A and B. 



the theoretical predictions and the simulation results. Moreover, the best fits obtained vary- 
ing the parameters A and B in Eq. (|T8|) are also plotted. Introducing density corrections to 
the expression for the temperature as discussed above for the pressure, is far from trivial. 
The arguments in Sec. [TTl leading to a closed separated equation for the temperature do not 
apply anymore and the resulting differential equations do not seem to have an elementary 
solution. In any case, the above results clearly confirm the accuracy of both the hydrody- 
namic equations and the used boundary conditions. Notice that the boundary layer next to 
the movable piston is much narrower for the temperature than for the pressure. 

Similar degree of agreement has been found for all the studied combinations of parameters 
in the intervals 20 < M/m < 120, 0.9 < a < 1, and 0.6 < «p < 1. In particular, the 
inelasticity of the collisions between the movable piston and the grains, measured by the 
coefficient ap, seems to affect very weakly the accuracy of the theory. 
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It is interesting to plot also the hydrodynamic profiles in the original scale z or, equiv- 
alently as a function of z* = zgo/v^. In Fig. [U the temperature profiles in the z*-scale 
are shown for the same two steady states as in the previous figures. Also included are the 
profiles obtained with the best fitting parameters A and B. It is observed that the agree- 
ment between theory and simulation is now worse than in the previous figures, where the 
spatial scale £ was employed. This is not at all surprising, since in the transformation of 
the theoretical prediction from the scale £ to the scale z, Eq. (fit)]) , the density profile is 
involved, and the ideal gas equation of state has been used. Consequently, the density effects 
discussed above manifest themselves in the transformation. Moreover, the scale z is defined 
from £ in a cumulative way, as an integral over the density profile, so that the discrepancies 
increase as z decreases. This effect is also clearly identified in Fig. where the profiles 
for the reduced density n* are given. It is also seen that the simulation results extend to 
larger values of z than the theoretical predictions for the profiles. It is worth to emphasize 
that the maximum value of the density, n max , remains quite low in both systems. For the 
system with a = 0.94, it is n max ~ 0.045a -2 ~ 0.039n c , while for the one with a = 0.98, 
n max ~ 0.04a- 2 ~ 0.035n c . 

The previous analysis of the profiles in the £ variable suggests that the analytical ex- 
pression of the reduced temperature in that variable is rather robust, in the sense of being 
very little affected by the finite density effects that, on the other hand, produce identifiable 
modifications on the pressure and density profiles. To check this idea, which can be useful 
to describe real experiments, the simulation profiles in Figs. H] and [5] have been also fitted as 
follows. The pressure is given by Eq. (TT5]) . and for the temperature, T(£), Eq. ( fl~8l) is used, 
with A and B being adjustable parameters. Nevertheless, the equation of state is not that 
of an ideal gas, but Eq. fl82l) . This equation of state is used both to compute the pressure 
from the simulation data, as discussed above, and also to transform from the £ coordinate to 
the z one by means of Eq. ( flOl) . The profiles obtained in this way are also included in Figs. 
H] and [5J and reproduce fairly well the hydrodynamic profiles in the bulk on the system. 

B. Global properties and boundary conditions 

The theoretical result for the average position of the piston L is given by Eq. flSUl) . with 
the density profile given by Eq. ( 1201) and the values of the constants A and B following from 
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FIG. 4: (Color online) Temperature profiles for the same systems as in Fig. [21 but now as a function 
of the reduced original coordinate z* defined in the text. The solid line is again the theoretical 
prediction with no adjustable parameters, and the symbols MD simulation results. The dotted- 
dashed lines (green online) have been obtained by fitting the two parameters A and B appearing 
in the expression of the temperature field using, moreover, the equation of state ([82]) . as discussed 
in the main text. 

Eqs. ( ITo]) and f|77|) . In Fig. [HI the above prediction is compared with MD simulation results. 
The value of the dimensionless height L* = Lgo/v^y as a function of M/m is plotted for two 
systems, one with a = 0.94, ap = 0.99, and the other one with a = 0.99, ap = 0.8. The 
error bars in the simulations data have been obtained from the mean square fluctuations of 
the position of the piston around its average value, once in the steady state. Although there 
is a systematic underestimation of L, the agreement between theory and simulation is good 
over a quite wide range of values of the mass ratio. Prompted by the previous discussion 
of the hydrodynamic profiles, it is tempting to investigate whether some of the observed 
discrepancy is due to finite density corrections. Then, instead of using Eq. (120]) for the 
density profile, the equation of state (182!) was employed. In the latter, the pressure was given 
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FIG. 5: (Color online) Density profiles for the same systems as in Fig. [2] as a function of the 
reduced coordinate z* defined in the text. The solid line is the theoretical prediction with no 
adjustable parameters and the symbols MD simulation results. The dotted-dashed lines (green 
online) have been obtained by fitting the two parameters A and B appearing in the expression of 
the temperature field using, moreover, the equation of state (f82j) . as discussed in the main text. 

by Eq. (I15p and the temperature by Eq. (1181) . with A and B still determined from Eqs. (175]) 
and (177j) . The problem is that now the resulting expression for L* depends on the velocity 
of the bottom vibrating wall, vw, and the simulation data in Fig. ([6]) have been obtained 
with different values of this velocity. For this reason, two different results are reported in the 
figure, corresponding to the two extreme values of v w used in the simulations. Although the 
dependence of the results on v\y is small, it is clearly appreciable on the scale of the figure. 
As expected, the curve corresponding to the largest value of v w is the closest to the ideal 
gas prediction, since the larger vw the more dilute the system. Although some discrepancy 
between theory and simulation still persists, increasing as the value of a decreases, it seems 
fair to conclude the presence of finite density effects in the simulation results. 
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FIG. 6: Average position of the piston L* , measured in the dimensionless units defined in the 
text, as a function of the the mass ratio M/m, for two different pairs of the restitution coefficients 
a, ap, as indicated in the figures. The symbols are MD simulation results and the solid lines the 
theoretical prediction derived in Sec. HVl The dashed lines and the dot-dashed lines have been 
obtained by using the equation of state (|82p . and correspond to the largest and smallest values of 
the velocity of the bottom wall vw used in the simulation, respectively. 

Another relevant property characterizing globally the system is the power dissipated per 
unit of section W, as a consequence of the inelasticity of collisions. As discussed at the end 
of Sec. HV] this quantity is related with the heat flux at the boundaries by 

D = Q L ~ Qo. (83) 

The fluxes Ql and Q have been computed in Sec. [TV] by using approximate kinetic theory 
arguments and the expressions derived for them are given by Eqs. (|66[) and (|68p . respectively. 
These expressions are closed by means of Eqs. (lOOj) and (1701) and making, as above, the 
additional approximation Tp/Tl — 1 in the expressions of 0, Eq. (167|) . 

In the MD simulations, D is directly obtained from the loss of kinetic energy in each 
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FIG. 7: Dimensionless dissipated power D* per unit of transversal length in the steady state, as 
a function of the ratio between the mass of the piston, M, and the mass of the grains, m. The 
symbols are MD simulation results, while the lines have been obtained from the expressions of the 
fluxes at the boundaries derived in this paper. The circles and solid line correspond to systems 
with a = 0.98 and ap = 0.8, and the squares and the dashed line to a = 0.94 and ap = 0.99. 

collision. In Fig. [7J 

D* = (84) 
mg v w 

is plotted as a function of M/m for the same pairs of values a, ap as in Fig. El The agreement 
between theory and MD simulation is quite good. This provides strong additional support 
for the analysis of the boundaries carried out in Sees. II I II and IIVI and, consequently, for the 
boundary conditions employed in the hydrodynamic equations. 



VI. CONCLUSIONS 



In this paper, the steady state of a fmidized granular gas of smooth inelastic hard spheres 
in presence of gravity, and with a movable lid on the top, has been investigated. It has been 
shown that, for weak inelasticity of the grain-grain collisions, the inelastic Navier-Stokes 
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equations provide an accurate description of the density, temperature, and pressure profiles. 
For stronger inelasticity of the grains, a hydrodynamic description beyond the Navier-Stokes 
approximation is needed due to the coupling between gradients and inelasticity. 

The issue of the required boundary conditions to solve the hydrodynamic equations has 
also been addressed. Here, the system was ffuidized by means of a vibrating elastic wall at the 
bottom. For simplicity, it was assumed to move in a sawtooth way and with high frequency 
and small amplitude. The collisions of the particles with the piston on the top were modeled 
as smooth and inelastic. The above mechanical characterization of the walls was translated 
into exact boundary conditions for the joint distribution functions of each of the walls and 
the gas next to it. These kinetic boundary conditions are general, with restriction neither 
in the gas density nor in the inelasticity. They must be taken into account when finding 
approximate hydrodynamic boundary conditions. Otherwise, relevant exact properties, as 
the hydrodynamic meaning of the pressure of the gas at the boundaries, can be violated. 
Here, an approximation consistent with the kinetic boundary condition was developed and 
shown to lead to quite accurate results. It must be emphasized that the aim here has not 
been to analyze in detail the boundary layers next to the walls, but rather to identify the 
effective boundary conditions that must be satisfied by the hydrodynamic fields extrapolated 
to the boundaries from the bulk of the system. 

Very recently, the so-called Knudsen temperature jump at a thermal wall has been in- 
corporated at the Navier-Stokes equations of a weakly inelastic dilute gas of hard disks [34]. 
As mentioned in the Introduction, thermal walls have not been proven to correspond to any 
well defined limit of vibrating walls. For the specific case of the inelastic piston on top of the 
system considered here, any analogy with a thermal wall seems very hard to justify. More- 



over, the approximation followed in 



34J involves the determination of a constant prefactor 



by means of MD simulations, while here all the constants are determined by the theory. 

The accuracy of the theoretical predictions has been checked by comparison with MD 
simulations. For some of the properties measured, e.g. the pressure profile, finite density 
effects have been identified and a practical way of accounting for them has been proposed. 
Of course, these effects could have been avoided by using the direct simulation Monte Carlo 
(DSMC) method [h], [jsj], instead of molecular dynamics. Nevertheless, the latter seems 
more appropriate in the present context, in which the interest focuses on the hydrodynamic 
description of a given state and not on a property of the Boltzmann equation. MD simula- 
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tions provide results closer to real experiments since they are not based on the validity of 
any kinetic equation. 

The interest has been on the description of the granular gas and its interaction with the 
piston. Of course, the study of the properties of the piston itself is also of great interest. 
For instance, its position fluctuations give the volume fluctuations of the gas. Some partial 
results have already been reported elsewhere [32j, but many points still deserve additional 
attention. 
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APPENDIX A: TRANSPORT COEFFICIENTS FOR A DILUTE GRANULAR 
GAS 



Here, the expressions of the dimensionless functions introduced in Sec. [TT] are given for 

nn 

the sake of completeness. They read [3, 

2d 
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These results have been obtained from the inelastic Boltzmann equation for hard spheres 
and disks by means of a generalized Chapman-Enskog algorithm, in the so-called first Sonine 
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approximation 
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APPENDIX B: CALCULATION OF THE FORCE EXERTED BY THE GAS ON 
THE PISTON IN THE STEADY STATE 



The force exerted by the granular gas on the piston equals the rate of momentum trans- 
ferred from the gas to the piston in the collisions. When a particle collides with the piston, 
the amount of momentum given by the former to the latter is given by Eq. (1261) . Then, the 
average force on the piston F z is 

/poo poo 
dx dZ dV z <d(g z )g z AP z $ 0:St (x, Z, V Z )5(Z - z) 
JO J -oo 



OO POO 



in M 

—— (1 + a P ) I dx I dZ I dV z g 2 z $ +ySt (x, Z, V Z )5(Z - z). (Bl) 
m + M 



J -oo 



This expression can be rewritten in a more familiar form. First, decompose it as 

//•oo poo 
dx I dZ dV z v z g z $ + , st (x,Z,V z )5(Z - z) 
JO J -oo 

f f°° f°° \ M 1 

-m dx dZ dV z v z — — (l + a P )g z g z $ + , st (x, Z, V Z )S(Z - z). 

J Jo J-oo L m + M J 

(B2) 

Next, taking into account Eq. (1231) . the second term on the right hand side of this equality 
is seen to be equivalent to 

/roo roo 
dx dZ dV z v' z g z $ +tgt (x, Z, V Z )5(Z - z) 
JO J - oc 

/POO POO 
dx* dZ dV* v z g*^ +iSt (r, v*, Z, V Z )5{Z - z) 
Jo J -oo 

OO POO 



= m dx dZ dV z v z g z $- >st (x,Z,V z )5(Z-z). (B3) 

J Jo J -oo 

In the above transformations, x* = {r,v*} and the exact boundary condition in Eq. (|4T|) 
has been employed. Substitution of Eq. flB3j) into Eq. (IB2I) and use of the parity of $o,st 
with respect to V z as well as its independence from r± gives 

/POO POO 
dx dZ dV z v z g z ^ 0jSt (x,Z,V z )S(Z - z) 
Jo J -oo 

= Wn L T L>z , (B4) 

where tll is the density of the granular gas next to the piston and T^ j2 a temperature 
parameter defined from the ^-component of the velocity. Its expression is given in Eq. ( |48l ). 

It is worth to stress that no assumption has been made over the properties of $o,st other 
than the associated with the general symmetry properties of the steady macroscopic state 
under consideration and the boundary condition fj4T|) . 
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APPENDIX C: THE BOUNDARY CONDITION AT THE VIBRATING WALL 



Accordingly with the description of the vibrating wall located at the bottom of the system 
given at the beginning of Sec. IIII[ when a particle with velocity v, being v z < 0, collides 
with this wall, its velocity is instantaneously changed into v' defined by 

v' z = 2v w - v z , v' ± = v ±} (CI) 

Here, as already mentioned, vyy is the (upwards) velocity of the wall and v± is the vector 
component of v perpendicular to the z-axis. The kinetic energy gained by the particle in 
such a collision is 

Ae = 2m(v^-v w v z ). (C2) 

For physical initial conditions in which there are no particles below the vibrating wall, the 
one-particle distributions function of the gas at arbitrary later times verifies 

f(x,t) = e(z)f (x,t), (C3) 

where fo(x, t) and its derivatives can be taken regular at z — 0. Next, the above distribution 
function is decomposed at the wall as 

f (x,t)6(z) = fj ) + \x,t)6(z) + ft ) (x 1 t)5(z), (C4) 

with 

ft\x,t) = Q(±v z )f (x,t). (C5) 
Conservation of the flux of particles at this wall is expressed by 

ft\x,t)\v z \S(z)dv = f!> +) (x',t)\v' z \5(z)dv', (C6) 

where x' = {r, v'}. Using Eqs. fIClj) . the above relation can be transformed into 

ri + \x,t)v z 6(z) = ft\x / ,t)\v>\6(z)- (C7) 

The heat flux at the vibrating wall, Qo, is given by the rate of energy input per unit of 
section of the wall, i.e. 

Qo{r±,t) = J dz J dv \v z \Aeft\x,t)5{z) 

= J dz J vmv 2 z v w fQ~\x,t)5(z)5(z) 

+ J dz J dvm\v z \vw\2vw — Vz\fo'\x,t)5(z). (C8) 
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Consider the second term on the right hand side of the above expression. Use of Eq. (1C7|) 
yields 

J ' dz J dvm\v z \vw\^vw — v z \fo~\x,t)5(z) = J ' dz J dv mv 2 z vw fo + \x,t)5(z), (C9) 
and substitution of this into Eq. (1C8P gives 

Qa{r±,t) = j dz J dvmv 2 z v w fo(x,t)5(z) = v w n (r ± ,t)T 0tZ (r ± ,t), (CIO) 



where 



and 



n (r ± ,t) = J dz J dvf (x,t)5(z) (Gil) 



n {r±,t)T 0>z {r ± ,t) = J dz J dvmv 2 J (x,t)5(z). (C12) 
Equation fl68|) in the main text follows by particularizing Eq. (IClOj) for the one-dimensional 



steady state considered in this paper and assuming isotropy of the pressure tensor of the gas 
next to the vibrating wall. 

The issue of the energy injected in a granular gas by means of an elastic vibrating wall 



moving in a sawtooth way has been addressed previously in ref. 28|. The arguments there 
are restricted to a particular state with the gas modeled by means of Gaussian distributions. 
Moreover, although the results reported in [28( are consistent with Eq. (IClOp . the simplicity 
of this latter form of expressing the result was not realized there. 
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